#ifdef __MKL__

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/storage.hpp>

#include "mklwrapper.h"

// #define __PRINT_STATS__

/*
********************************************************************************
*                              INTEL CONFIDENTIAL
*   Copyright(C) 2004-2010 Intel Corporation. All Rights Reserved.
*   The source code contained  or  described herein and all documents related to
*   the source code ("Material") are owned by Intel Corporation or its suppliers
*   or licensors.  Title to the  Material remains with  Intel Corporation or its
*   suppliers and licensors. The Material contains trade secrets and proprietary
*   and  confidential  information of  Intel or its suppliers and licensors. The
*   Material  is  protected  by  worldwide  copyright  and trade secret laws and
*   treaty  provisions. No part of the Material may be used, copied, reproduced,
*   modified, published, uploaded, posted, transmitted, distributed or disclosed
*   in any way without Intel's prior express written permission.
*   No license  under any  patent, copyright, trade secret or other intellectual
*   property right is granted to or conferred upon you by disclosure or delivery
*   of the Materials,  either expressly, by implication, inducement, estoppel or
*   otherwise.  Any  license  under  such  intellectual property  rights must be
*   express and approved by Intel in writing.
*
********************************************************************************
*   Content : MKL PARDISO C example
*
********************************************************************************
*/
/* -------------------------------------------------------------------- */
/* Example program to show the use of the "PARDISO" routine */
/* on symmetric linear systems */
/* -------------------------------------------------------------------- */
/* This program can be downloaded from the following site: */
/* www.pardiso-project.org */
/* */
/* (C) Olaf Schenk, Department of Computer Science, */
/* University of Basel, Switzerland. */
/* Email: olaf.schenk@unibas.ch */
/* -------------------------------------------------------------------- */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <vector>
#include <iostream>

#include "mkl_pardiso.h"
#include "mkl_types.h"

// #define __PRINT_STATS__

void SolvePardiso(const boost::numeric::ublas::mapped_matrix<double>& A, const boost::numeric::ublas::vector<double, std::vector<double> >& vb,  boost::numeric::ublas::vector<double, std::vector<double> >& vx, const bool symmetric) {
	namespace ublas = boost::numeric::ublas;
	typedef const ublas::mapped_matrix<double> SparseMatrix;
	typedef ublas::matrix_row< SparseMatrix > Row;
	
	std::vector<MKL_INT> via, vja;
	std::vector<double> va;
	MKL_INT k = 0;
	for(size_t i = 0; i < A.size1(); i++) {
		via.push_back(k + 1);
		Row r(A, i);
		for(Row::iterator itr = r.begin(); itr != r.end(); itr++) {
			const size_t j = itr.index();
			const double d = *itr;
			if((symmetric && j >= i) || !symmetric) {
				va.push_back(d);
				vja.push_back(j + 1);
				k++;
			}
		}
	}
	via.push_back(k + 1);
	
#ifdef __PRINT_STATS__
// 	std::cout << "ia: "; for(size_t i = 0; i < via.size(); i++) std::cout << via[i] << " "; std::cout << std::endl;
// 	std::cout << "ja: "; for(size_t i = 0; i < vja.size(); i++) std::cout << vja[i] << " "; std::cout << std::endl;
	std::cout << "a: "; for(size_t i = 0; i < va.size(); i++) {std::cout << "va: " << va[i] << " "; std::cout << std::endl;}
#endif

	MKL_INT* ia = &via[0];
	MKL_INT* ja = &vja[0];
	double* a = &va[0];
	
	MKL_INT n = A.size1();

	/* Matrix data. */
	/*MKL_INT n = 8;
	MKL_INT ia[ 9] = { 1, 5, 8, 10, 12, 15, 17, 18, 19 };
	MKL_INT ja[18] = { 1, 3, 6, 7,
		2, 3, 5,
		3, 8,
		4, 7,
		5, 6, 7,
		6, 8,
		7,
		8 };
	double a[18] = { 7.0, 1.0, 2.0, 7.0,
		-4.0, 8.0, 2.0,
		1.0, 5.0,
		7.0, 9.0,
		5.0, 1.0, 5.0,
		-1.0, 5.0,
		11.0,
		5.0 };*/

	MKL_INT mtype;
	if(symmetric)
		mtype = -2; /* Real symmetric matrix */
	else
		mtype = 11; /* Real and nonsymmetric matrix */
	/* RHS and solution vectors. */
	std::vector<double> vx2(n);
	double* x = &vx2[0];
	MKL_INT nrhs = 1; /* Number of right hand sides. */
	/* Internal solver memory pointer pt, */
	/* 32-bit: int pt[64]; 64-bit: long int pt[64] */
	/* or void *pt[64] should be OK on both architectures */
	void *pt[64];
	/* Pardiso control parameters. */
	MKL_INT iparm[64];
	MKL_INT maxfct, mnum, phase, error, msglvl;
	/* Auxiliary variables. */
	MKL_INT i;
	double ddum; /* Double dummy */
	MKL_INT idum; /* Integer dummy. */
/* -------------------------------------------------------------------- */
/* .. Setup Pardiso control parameters. */
/* -------------------------------------------------------------------- */
	for (i = 0; i < 64; i++) {
		iparm[i] = 0;
	}
	iparm[0] = 1; /* No solver default */
	iparm[1] = 2; /* Fill-in reordering from METIS */
	/* Numbers of processors, value of OMP_NUM_THREADS */
	iparm[2] = 1;
	iparm[3] = 0; /* No iterative-direct algorithm */
	iparm[4] = 0; /* No user fill-in reducing permutation */
	iparm[5] = 0; /* Write solution into x */
	iparm[6] = 0; /* Not in use */
	iparm[7] = 2; /* Max numbers of iterative refinement steps */
	iparm[8] = 0; /* Not in use */
	iparm[9] = 13; /* Perturb the pivot elements with 1E-13 */
	iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */
	iparm[11] = 0; /* Not in use */
	iparm[12] = 0; /* Maximum weighted matching algorithm is switched-off (default for symmetric). Try iparm[12] = 1 in case of inappropriate accuracy */
	iparm[13] = 0; /* Output: Number of perturbed pivots */
	iparm[14] = 0; /* Not in use */
	iparm[15] = 0; /* Not in use */
	iparm[16] = 0; /* Not in use */
	iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
	iparm[18] = -1; /* Output: Mflops for LU factorization */
	iparm[19] = 0; /* Output: Numbers of CG Iterations */
	maxfct = 1; /* Maximum number of numerical factorizations. */
	mnum = 1; /* Which factorization to use. */

#ifdef __PRINT_STATS__
	msglvl = 1; /* Print statistical information in file */
#else
	msglvl = 0; /* Print statistical information in file */
#endif
	error = 0; /* Initialize error flag */
/* -------------------------------------------------------------------- */
/* .. Initialize the internal solver memory pointer. This is only */
/* necessary for the FIRST call of the PARDISO solver. */
/* -------------------------------------------------------------------- */
	for (i = 0; i < 64; i++) {
		pt[i] = 0;
	}
/* -------------------------------------------------------------------- */
/* .. Reordering and Symbolic Factorization. This step also allocates */
/* all memory that is necessary for the factorization. */
/* -------------------------------------------------------------------- */
	phase = 11;
	PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
		&n, a, ia, ja, &idum, &nrhs,
		iparm, &msglvl, &ddum, &ddum, &error);
	if (error != 0) {
		printf("\nERROR during symbolic factorization: %d", error);
		exit(1);
	}
#ifdef __PRINT_STATS__
	printf("\nReordering completed ... ");
	printf("\nNumber of nonzeros in factors = %d", iparm[17]);
	printf("\nNumber of factorization MFLOPS = %d", iparm[18]);
#endif
/* -------------------------------------------------------------------- */
/* .. Numerical factorization. */
/* -------------------------------------------------------------------- */
	phase = 22;
	PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
		&n, a, ia, ja, &idum, &nrhs,
		iparm, &msglvl, &ddum, &ddum, &error);
	if (error != 0) {
		printf("\nERROR during numerical factorization: %d", error);
		exit(2);
	}
#ifdef __PRINT_STATS__
	printf("\nFactorization completed ... ");
#endif
/* -------------------------------------------------------------------- */
/* .. Back substitution and iterative refinement. */
/* -------------------------------------------------------------------- */
	phase = 33;
	iparm[7] = 2; /* Max numbers of iterative refinement steps. */
	
	std::vector<double> vb2(n);
	for(size_t i = 0; i < n; i++)
		vb2[i] = vb(i);
	/* Set right hand side to one. */
// 	for (i = 0; i < n; i++) {
// 		b[i] = 1;
// 	}

	double* b = &vb2[0];
	PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
		&n, a, ia, ja, &idum, &nrhs,
		iparm, &msglvl, b, x, &error);
	if (error != 0) {
		printf("\nERROR during solution: %d", error);
		exit(3);
	}
#ifdef __PRINT_STATS__
	printf("\nSolve completed ... ");
	printf("\nThe solution of the system is: ");
	for (i = 0; i < n; i++) {
		printf("\n x [%d] = % f", i, x[i] );
	}
	printf ("\n");
#endif

	for(size_t i = 0; i < n; i++) {
		vx(i) = vx2[i];
	}
/* -------------------------------------------------------------------- */
/* .. Termination and release of memory. */
/* -------------------------------------------------------------------- */
	phase = -1; /* Release internal memory. */
	PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
		&n, &ddum, ia, ja, &idum, &nrhs,
		iparm, &msglvl, &ddum, &ddum, &error);
}

void SolvePardiso(const boost::numeric::ublas::mapped_matrix<double>& A, const std::vector< boost::numeric::ublas::vector<double, std::vector<double> > >& vb,  std::vector< boost::numeric::ublas::vector<double, std::vector<double> > >& vx, const bool symmetric) {
	namespace ublas = boost::numeric::ublas;
	typedef const ublas::mapped_matrix<double> SparseMatrix;
	typedef ublas::matrix_row< SparseMatrix > Row;
	
	std::vector<MKL_INT> via, vja;
	std::vector<double> va;
	MKL_INT k = 0;
	for(size_t i = 0; i < A.size1(); i++) {
		via.push_back(k + 1);
		Row r(A, i);
		for(Row::iterator itr = r.begin(); itr != r.end(); itr++) {
			const size_t j = itr.index();
			const double d = *itr;
			if((symmetric && j >= i) || !symmetric) {
				va.push_back(d);
				vja.push_back(j + 1);
				k++;
			}
		}
	}
	via.push_back(k + 1);
	
#ifdef __PRINT_STATS__
	std::cout << "ia: "; for(size_t i = 0; i < via.size(); i++) std::cout << via[i] << " "; std::cout << std::endl;
	std::cout << "ja: "; for(size_t i = 0; i < vja.size(); i++) std::cout << vja[i] << " "; std::cout << std::endl;
	std::cout << "a: "; for(size_t i = 0; i < va.size(); i++) std::cout << va[i] << " "; std::cout << std::endl;
#endif

	MKL_INT* ia = &via[0];
	MKL_INT* ja = &vja[0];
	double* a = &va[0];
	
	MKL_INT n = A.size1();

	/* Matrix data. */
	/*MKL_INT n = 8;
	MKL_INT ia[ 9] = { 1, 5, 8, 10, 12, 15, 17, 18, 19 };
	MKL_INT ja[18] = { 1, 3, 6, 7,
		2, 3, 5,
		3, 8,
		4, 7,
		5, 6, 7,
		6, 8,
		7,
		8 };
	double a[18] = { 7.0, 1.0, 2.0, 7.0,
		-4.0, 8.0, 2.0,
		1.0, 5.0,
		7.0, 9.0,
		5.0, 1.0, 5.0,
		-1.0, 5.0,
		11.0,
		5.0 };*/

	MKL_INT mtype;
	if(symmetric)
		mtype = -2; /* Real symmetric matrix */
	else
		mtype = 11; /* Real and nonsymmetric matrix */
	/* RHS and solution vectors. */
	std::vector<double> vx2(n);
	double* x = &vx2[0];
	MKL_INT nrhs = 1; /* Number of right hand sides. */
	/* Internal solver memory pointer pt, */
	/* 32-bit: int pt[64]; 64-bit: long int pt[64] */
	/* or void *pt[64] should be OK on both architectures */
	void *pt[64];
	/* Pardiso control parameters. */
	MKL_INT iparm[64];
	MKL_INT maxfct, mnum, phase, error, msglvl;
	/* Auxiliary variables. */
	MKL_INT i;
	double ddum; /* Double dummy */
	MKL_INT idum; /* Integer dummy. */
/* -------------------------------------------------------------------- */
/* .. Setup Pardiso control parameters. */
/* -------------------------------------------------------------------- */
	for (i = 0; i < 64; i++) {
		iparm[i] = 0;
	}
	iparm[0] = 1; /* No solver default */
	iparm[1] = 2; /* Fill-in reordering from METIS */
	/* Numbers of processors, value of OMP_NUM_THREADS */
	iparm[2] = 1;
	iparm[3] = 0; /* No iterative-direct algorithm */
	iparm[4] = 0; /* No user fill-in reducing permutation */
	iparm[5] = 0; /* Write solution into x */
	iparm[6] = 0; /* Not in use */
	iparm[7] = 2; /* Max numbers of iterative refinement steps */
	iparm[8] = 0; /* Not in use */
	iparm[9] = 13; /* Perturb the pivot elements with 1E-13 */
	iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */
	iparm[11] = 0; /* Not in use */
	iparm[12] = 0; /* Maximum weighted matching algorithm is switched-off (default for symmetric). Try iparm[12] = 1 in case of inappropriate accuracy */
	iparm[13] = 0; /* Output: Number of perturbed pivots */
	iparm[14] = 0; /* Not in use */
	iparm[15] = 0; /* Not in use */
	iparm[16] = 0; /* Not in use */
	iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
	iparm[18] = -1; /* Output: Mflops for LU factorization */
	iparm[19] = 0; /* Output: Numbers of CG Iterations */
	maxfct = 1; /* Maximum number of numerical factorizations. */
	mnum = 1; /* Which factorization to use. */

#ifdef __PRINT_STATS__
	msglvl = 1; /* Print statistical information in file */
#else
	msglvl = 0; /* Print statistical information in file */
#endif
	error = 0; /* Initialize error flag */
/* -------------------------------------------------------------------- */
/* .. Initialize the internal solver memory pointer. This is only */
/* necessary for the FIRST call of the PARDISO solver. */
/* -------------------------------------------------------------------- */
	for (i = 0; i < 64; i++) {
		pt[i] = 0;
	}
/* -------------------------------------------------------------------- */
/* .. Reordering and Symbolic Factorization. This step also allocates */
/* all memory that is necessary for the factorization. */
/* -------------------------------------------------------------------- */
	phase = 11;
	PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
		&n, a, ia, ja, &idum, &nrhs,
		iparm, &msglvl, &ddum, &ddum, &error);
	if (error != 0) {
		printf("\nERROR during symbolic factorization: %d", error);
		exit(1);
	}
#ifdef __PRINT_STATS__
	printf("\nReordering completed ... ");
	printf("\nNumber of nonzeros in factors = %d", iparm[17]);
	printf("\nNumber of factorization MFLOPS = %d", iparm[18]);
#endif
/* -------------------------------------------------------------------- */
/* .. Numerical factorization. */
/* -------------------------------------------------------------------- */
	phase = 22;
	PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
		&n, a, ia, ja, &idum, &nrhs,
		iparm, &msglvl, &ddum, &ddum, &error);
	if (error != 0) {
		printf("\nERROR during numerical factorization: %d", error);
		exit(2);
	}
#ifdef __PRINT_STATS__
	printf("\nFactorization completed ... ");
#endif
/* -------------------------------------------------------------------- */
/* .. Back substitution and iterative refinement. */
/* -------------------------------------------------------------------- */
	phase = 33;
	iparm[7] = 2; /* Max numbers of iterative refinement steps. */
	
	for(size_t k = 0; k < vb.size(); k++) {
		std::vector<double> vb2(n);
		for(size_t i = 0; i < n; i++)
			vb2[i] = vb[k](i);
		/* Set right hand side to one. */
		// 	for (i = 0; i < n; i++) {
		// 		b[i] = 1;
		// 	}

		double* b = &vb2[0];
		PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
				&n, a, ia, ja, &idum, &nrhs,
				iparm, &msglvl, b, x, &error);
		if (error != 0) {
			printf("\nERROR during solution: %d", error);
			exit(3);
		}
#ifdef __PRINT_STATS__
		printf("\nSolve completed ... ");
		printf("\nThe solution of the system is: ");
		for (i = 0; i < n; i++) {
			printf("\n x [%d] = % f", i, x[i] );
		}
		printf ("\n");
#endif

		for(size_t i = 0; i < n; i++) {
			vx[k](i) = vx2[i];
		}
	}
/* -------------------------------------------------------------------- */
/* .. Termination and release of memory. */
/* -------------------------------------------------------------------- */
	phase = -1; /* Release internal memory. */
	PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
		&n, &ddum, ia, ja, &idum, &nrhs,
		iparm, &msglvl, &ddum, &ddum, &error);
}
#endif
